$onText
DISE-2024 (Dynamic Integrated Space-Economy) model
File: DISE_2024_04.gms (Version 0.3)
Decreasing launch costs
August 2025
Reference: Optimal path for orbital debris.
Bongers, A. and Torres, J. L. (2025).
$offText


*	Simulation period
Set t	'Time periods'				/2023*2272/;

Sets
        tfirst(t)   'First period (2023)'
        tlast(t)    'Last period (T)'
    	tnotlast(t) 'All periods but last'
	    ;
        
        tfirst(t)$(ord(t)=1)       =   yes;
        tlast(t)$(ord(t)=card(t))  =   yes;
        tnotlast(t)                =   not tlast(t);

Parameters

* 	Economic parameters

** 	Preferences parameters
	rho	'Intertemporal pure preference rate'		/0.015/
	sigma	'Relative risk-aversion'			/1.5/
    
**  Exogenous growth sources
    	gN0	    'Labor growth rate parameter'			/0.05/
    	Nstar	'Asymptotic population'				/10200/
    	gA0	    'TFP growth rate'				        /0.015/
    	deltaA	'TFP growth rate trend'				/0.001/
    	gQ0	    'Satellite ISTC growth rate'			/0.03/
    	deltaQ	'Satellite ISTC growth rate trend'	/0.005/

**	Initial values economic variables
	    k0	    'Initial Earth capital'			/555.6987/         		
    	s0     	'Initial Satellite capital'			/1.1959/  			
    	y0	    'Initial output'				/184.65/
    	x0      'Initial value of satellites destroyed'
    	N0      'Initial population'				/8056/
        
** 	Technological parameters      
	    alpha1 	'Output-capital elasticity'			/0.3479/
    	alpha2 	'Output-satellite elasticity'			/0.0021/ 
    	deltak 	'Capital depreciation rate'			/0.07/
    	deltas 	'Satellite depreciation rate'			/0.15/ 
	    a0	    'Initial value of TFP'				      
	    q0      'Initial value of satellite ISTC'		/1.00/

**  Launch cost
	b0      'Initial proportion of launch cost'     /0.3/
    gb0     'Initial change in launch cost'         /-0.05/
    deltab  'Launch cost growth rate trend'         /0.1/

** 	Physical parameters
	deltaf	'Fragments decay rate'				/0.01/
	deltaw  'Dereclit satellites decay rate'        	/0.00015/
	deltaz  'Rocket bodies decay rate'			/0.00015/
	omega 	'Number of fragments per launch'		/4.00/
	xi	    'Fraction of derelict satellites in orbit'	/0.40/
	theta	'Collision risk'				/1.25e-10/
	epsilonw'Fraction of dead satellite breakups'		/0.001/
	epsilonz'Fraction of body rockets breakups'		/0.0012/
	psi	    'Rocket bodies per launch'			/0.60/
	phiw 	'Fragments per derelict satellite breakup'	/44.60/
	phiz    'Fragments per rocket body breakup'		/100.20/
	gammas	'Fragments per operational satellite collision'	/70/
	gammaw  'Fragments per derelict satellite collision'	/70/
	gammaz  'Fragments per rocket body collision'		/70/
	eta	    'Number of satellites per launch'		/13.6/
	v	    'Collision avoidance'				/0/
	mu	    'Mapping parameter'
	
**	Initial values physical variables	
	w0	    'Abandoned derelict satellites 2023'		/3500/
	z0	    'Rocket bodies 2023'				/2050/
	f0	    'Number of fragments larger than 10cm 2023'	/30000/
	d0	    'Total number of debris larger than 10cm 2023'	
	ss0	    'Operational satellites 2023'			/8500/
	xx0	    'Satellites destroyed by collisions 2023'
	ff0	    'Total number of fragments larger than 1cm 2023'
	fff0 	'Total number of fragments larger than 1mm 2023'
	dd0	    'Total number of debris larger than 1cm 2023'
	ddd0 	'Total number of debris larger than 1mm 2023'
	ob1	    'Number of objects larger than 10cm (ESA)'	/30000/
	ob2	    'Number of objects larger than 1cm (ESA)'	/1000000/
	ob3	    'Number of objects larger than 1mm (ESA)'	/130000000/
	nu1	    'Ratio >1cm to >10cm'
	nu2	    'Ratio >1mm to >1cm'
	nu3	    'Ratio >1mm to >10cm'
	;

*	Calculus of key constants and initial values
    	mu	=	ss0/s0;
    	nu1 	=   	ob2/ob1;
    	nu2 	=   	ob3/ob2;
    	nu3 	=   	ob3/ob1;
    	d0 	= 	w0+z0+f0;
	ff0	=	(1+nu1)*f0;
	dd0	= 	w0+z0+ff0;
	fff0	=	nu3*f0;
	ddd0	=	w0+z0+fff0;
	a0	=	y0/(k0**alpha1*s0**alpha2*N0**(1-alpha1-alpha2));
	x0	=	theta*dd0*s0;
	xx0 	=   	theta*dd0*ss0;
	

*   	Timing

Parameters
    	beta(t)	'Discount factor'
    	tval(t)	'numerical value of t';
	    tval(t)	=	ord(t)-1;

Parameters
*   Technology
        A(t)    'TFP'
        Q(t)    'Satellite ISTC'
        gA(t)   'Growth rate of TFG'
        gQ(t)   'Growth rate of Satellite ISTC'
        ;
    
        gA(t)       =   gA0*exp(-deltaA*(tval(t)-1));
        gQ(t)       =   gQ0*exp(-deltaQ*(tval(t)-1));
        a("2023")       =   a0;
        q("2023")       =   q0;
        loop(t,A(t+1)   =   a(t)*exp(gA(t)););
        loop(t,Q(t+1)   =   q(t)*exp(gQ(t)););
Parameters
        b(t)    'Fraction of launch cost'
        gb(t)   'Decline in the launch cost';
        
        gb(t)           =   gb0*exp(-deltab*(tval(t)-1));
        b("2023")       =   b0;
        loop(t,b(t+1)   =   b(t)*exp(gb(t)););
        
    
*   Population
Parameters
    N(t)    'Population';
    N("2023")       =   N0;
    loop(t,N(t+1)   =   N(t)*(Nstar/N(t))**gN0;);

	
* The terminal weight beta(tlast) computation.
	beta(tnotlast(t))	=	power(1+rho,-tval(t));
	beta(tlast(t))		=	(1/rho)*power(1+rho,1-tval(t)-1);
	
Variables

* Economic variables
	SU      'Total summatory discounted utility'
	y(t)	'Output'
   	c(t)    'Consumption'
   	ratcy(t)'Ratio consumption-output'
   	ik(t)	'Investment in Earth capital'
   	is(t)	'Investment in Space capital'
   	h(t)    'Investment in satellite'
   	l(t)    'Launch cost'
   	k(t)    'Earch capital stock'
	s(t)	'Stock of satellites'
	x(t)	'Satellites destroyed by collisions'
	ii(t)   'Investment'
	hh(t)	'Number of new satellites inserted into orbit'
	U(t)	'Instantaneous utility'
	DU(t)   'Discounted utility'
	srk(t)	'Investment rate in Earth capital'
	srs(t)	'Investment rate in Space capital'

* Physical variables
	xx(t)	'Number of destroyed satellites'
	ss(t)   'Number of operational satellites (stock)'
	d(t)	'Number of debris larger than 10cm'
	dd(t)	'Number of debris larger than 1cm'
	ddd(t)	'Number of debris larger than 1mm'
	ll(t)	'Number of launches'
	ds(t)	'Value of derelict satellites per period'
	dss(t)	'Number of derelict satellites per period'
	dso(t)	'Number of derelict satellites abandoned in orbit per period'
	dsb(t)	'Number of derelict satellites breakups'
	dsd(t)	'Number of derelict satellites destroyed by collision'
	w(t)	'Number of derelict satellites abandoned in orbit (stock)'
	rb(t)	'Number of rocket bodies as a function of launches'
	rbb(t)	'Number of rocket bodies breakups'
	rbd(t)	'Number of rocket bodies destroyed by collision'
	z(t)	'Number of rocket bodies (stock)'
	f(t)	'Number of fragments larger than 10cm'
	ff(t)	'Number of fragments larger than 1cm'
	fff(t)	'Number of fragments larger than 1mm'
	;

Equations

* 	Economic variables
	utility(t)      'Instantaneous utility'
	sumdutility	'Summatory of discounted utility'
	production(t)   'Cobb-Douglas production function'
	allocation(t)   'Household choose between consumption and investments'
	totalinvest(t)	'Total investment'
	launch(t)       'Launch cost'
	satel(t)        'Satellite cost'
	accumulak(t) 	'Capital accumulation'
	accumulas(t)	'Satellite accumulation'
	ifinal(t)       'Minimal capital investment in final period'
	hfinal(t)	'Minimal satellite investment in final period'  
	iifinal(t)	'Minimal total investment in final period'
	destroyed(t)	'Value of satellites destroyed by collision'
	ratiocy(t)      'Consumption-output ratio'
	ratioiy(t)	'Investment rate in Earth capital'
	ratiohy(t)	'Investment rate in Space capital'

*	Physical variables
	satellites(t)	'Number of operational satellites'
	deadsat(t)	'Derelict satellites'
	deadsaton(t)	'Derelict satellites abandoned in orbit'
	accumuladso(t)	'Stock of derelict satellites in orbit'
	rocketbodies(t)	'Number of new rocket bodies in orbit'
	accumularb(t)	'Stock of rocket bodies in orbit'
	fragments(t)	'Stock of fragments larger than 10cm in orbit'
	fragments1(t)	'Stock of fragments larger than 1cm in orbit'
	fragments2(t)	'Stock of fragments larger than 1mm in orbit'
	debris(t)	'Stock of debris larger than 10cm'
	debris1(t)	'Stock of debris larger than 1cm'
	debris2(t)	'Stock of debris larger than 1mm'
	newsat(t)	'Number of new satellites inserted into orbit'
	nlaunches(t)	'Number of launches'
	ndeadsat(t)	'Number of derelict satellites per period'
	ndestroyed(t)	'Number of satellites destroyed by collision'
	deadsatb(t)	'Number of derelict satellites breakups'
	deadsatd(t)	'Number of derelict satellites destroyed by collision'
	rbbreakups(t)	'Number of rocket bodies breakups'
	rbdestroyed(t)	'Number of rocket bodies destroyed by collision'
	limit(t)           'Restriction to limit satellites destroyed'
	;


* Equations of the model

** Utility

	utility(t)..	U(t)		=e=	((c(t)/N(t))**(1-sigma)-1)/(1-sigma);
	sumdutility..	SU  	 	=e=	sum(t, beta(t)*U(t)*N(t));

** Production, consumption and investment
	production(t)..    		y(t) 		=e= 	a(t)*(k(t)**alpha1)*(s(t)**alpha2)*(N(t)**(1-alpha1-alpha2));
	allocation(t)..    		c(t) 		=e= 	y(t)-ik(t)-is(t);
	totalinvest(t)..		ii(t)		=e=	    ik(t)+is(t);
	launch(t)..             l(t)        =e=     b(t)*is(t);
	satel(t)..              h(t)        =e=     is(t)-l(t);
	accumulak(tnotlast(t))..  	k(t+1) 		=e=	(1-deltak)*k(t)+ik(t);
	accumulas(tnotlast(t))..	s(t+1)		=e=	(1-deltas)*s(t)+q(t)*(1-b(t))*is(t)-x(t);
	ifinal(tlast)..     		ik(tlast) 	=g= 	(0.02+deltak)*k(tlast);
	hfinal(tlast)..			is(tlast)	=g=	(0.02+deltas)*s(tlast)/q(tlast);
	iifinal(tlast)..		ii(tlast)	=e=	ik(tlast)+is(tlast);
	ratiocy(t)..			ratcy(t)	=e=	c(t)/y(t);
	ratioiy(t)..			srk(t)		=e=	ik(t)/y(t);
	ratiohy(t)..			srs(t)		=e=	is(t)/y(t);

** Damages
	destroyed(t)..			x(t)		=e=	(1-v)*theta*dd(t)*s(t);
	ndestroyed(t)..			xx(t)		=e=	(1-v)*theta*dd(t)*ss(t);

** Satellite sector
	nlaunches(t)..			ll(t)		=e=	mu*h(t)/eta;
	deadsat(t)..			ds(t)		=e=	deltas*s(t);
	ndeadsat(t)..			dss(t)		=e=	deltas*ss(t);
	newsat(t)..			    hh(t)		=e=	eta*ll(t);
	satellites(t)..			ss(t)		=e=	mu*s(t)/q(t);

** Objects in orbit
    	deadsaton(t)..			dso(t)		=e=	xi*deltas*ss(t);
    	accumuladso(tnotlast(t))..	w(t+1)  	=e= 	(1-deltaw)*w(t)-epsilonw*w(t)-theta*(dd(t)+ss(t))*w(t)+dso(t);
    	rocketbodies(t)..		rb(t)       	=e= 	psi*ll(t);
    	accumularb(tnotlast(t))..   	z(t+1)		=e=	(1-deltaz)*z(t)-epsilonz*z(t)-theta*(dd(t)+ss(t))*z(t)+rb(t);
    	fragments(tnotlast(t))..    	f(t+1)   	=e= 	(1-deltaf)*f(t)+omega*ll(t)+phiw*dsb(t)+phiz*rbb(t)+gammas*xx(t)+gammaw*theta*dd(t)*w(t)+gammaz*theta*dd(t)*z(t);
    	fragments1(t)..             	ff(t)    	=e= 	(1+nu1)*f(t);
    	fragments2(t)..             	fff(t)    	=e= 	(1+nu3)*f(t);
    	debris(t)..                 	d(t)		=e= 	w(t)+z(t)+f(t);
    	debris1(t)..                	dd(t)       	=e= 	w(t)+z(t)+ff(t);
    	debris2(t)..                	ddd(t)       	=e= 	w(t)+z(t)+fff(t);
	deadsatb(t)..			dsb(t)		=e=	epsilonw*w(t);
	deadsatd(t)..			dsd(t)		=e=	theta*(dd(t)+ss(t))*w(t);
	rbbreakups(t)..			rbb(t)		=e=	epsilonz*z(t);
	rbdestroyed(t)..		rbd(t)		=e=	theta*(dd(t)+ss(t))*z(t);
    limit(t)..                  x(t)        =l= s(t);

* Bounds
	k.lo(t) = 0.001;
	s.lo(t)	= 0.001;
	c.lo(t) = 0.001;
	y.lo(t) = 0.001;
	is.lo(t)= 0;
	w.lo(t) = 0;
	z.lo(t) = 0;

* Initial conditions
	k.fx(tfirst) 	= 	k0;
	s.fx(tfirst) 	= 	s0;
	y.fx(tfirst)    =   y0;
	d.fx(tfirst) 	= 	d0;
	w.fx(tfirst) 	= 	w0;
	z.fx(tfirst) 	= 	z0;
	f.fx(tfirst)	=	f0;

* Solution options
   	option iterlim = 99900;
   	option reslim  = 99999;
   	option solprint = on;
	option limrow = 0;
	option limcol = 0;
    
model DISE / all /;

solve DISE maximizing SU using nlp;
solve DISE maximizing SU using nlp;


* Export to Excel using GDX utilities
execute_unload "soldise01b.gdx"
execute '=gdx2xls soldise01b.gdx';

* Reusable launch vehicles
    psi = 0;

model DISE2 / all /;
solve DISE2 maximizing SU using nlp;
solve DISE2 maximizing SU using nlp;
solve DISE2 maximizing SU using nlp;
execute_unload "soldise02b.gdx"
execute '=gdx2xls soldise02b.gdx';
   
* De-orbiting
    xi = 0;
    
model DISE3 / all /;
solve DISE3 maximizing SU using nlp;
solve DISE3 maximizing SU using nlp;
execute_unload "soldise03b.gdx"
execute '=gdx2xls soldise03b.gdx';

* No-breakups
    psi = 0.60;
    xi = 0.40;
    epsilonw=0;
    epsilonz=0;

model DISE4 / all /;
solve DISE4 maximizing SU using nlp;
solve DISE4 maximizing SU using nlp;
execute_unload "soldise04b.gdx"
execute '=gdx2xls soldise04b.gdx';

* Debris-free launch systems
    epsilonw=0.001;
    epsilonz=0.0012;
    omega = 0;
    psi = 0;
    
model DISE5 / all /;
solve DISE5 maximizing SU using nlp;
solve DISE5 maximizing SU using nlp;
execute_unload "soldise05b.gdx"
execute '=gdx2xls soldise05b.gdx';

* Combination of De-orbiting+No-breakups+Debris-free
    xi=0;
    epsilonw = 0;
    epsilonz = 0;
    
model DISE6 / all /;
solve DISE6 maximizing SU using nlp;
solve DISE6 maximizing SU using nlp;
execute_unload "soldise06b.gdx"
execute '=gdx2xls soldise06b.gdx';

* No-collision
    epsilonw = 0.001;
    epsilonz = 0.0012;
    psi = 0.60;
    xi = 0.40;
    omega = 4;
    v = 1;
    
model DISE7 / all /;
solve DISE7 maximizing SU using nlp;
solve DISE7 maximizing SU using nlp;
execute_unload "soldise07b.gdx"
execute '=gdx2xls soldise07b.gdx';

* ESA zero debris
    psi = 0;
    xi = 0;
    omega =0;
    v = 1;
    
model DISE8 / all /;
solve DISE8 maximizing SU using nlp;
solve DISE8 maximizing SU using nlp;
execute_unload "soldise08b.gdx"
execute '=gdx2xls soldise08b.gdx';
